Quantifying microbubble clustering in turbulent flow from single-point measurements 



Enrico Calzavarini*,^ Thomas H. van den Berg/ Federico Toschi,^ and Detlef Lohse^"'^ 

^Department of Applied Physics, JMBC Burgers Center for Fluid Dynamics, 
and IMPACT Institute, University of Twente, 7500 AE Enschede, The Netherlands. 
^ lAC-CNR, Istituto per le Applicazioni del Calcolo, Viale del Policlinico 137, 
1-00161 Roma, Italy and INFN, via Saragat 1, 1-44100 Ferrara, Italy. 
(Dated: February 2, 2008) 

Single-point hot-wire measurements in the bulk of a turbulent channel have been performed in 
order to detect and quantify the phenomenon of preferential bubble accumulation. We show that 
statistical analysis of the bubble-probe coUiding-times series can give a robust method for investi- 
gation of clustering in the bulk regions of a turbulent flow where, due to the opacity of the flow, 
no imaging technique can be employed. We demonstrate that microbubbles {Ro — 100 fim) in a 
developed turbulent flow, where the Kolmogorov length-scale is 77 ~ _Ro, display preferential con- 
centration in small scale structures with a typical statistical signature ranging from the dissipative 
range, 0{rj), up to the low inertial range O{100r]). A comparison with Eulerian-Lagrangian numeri- 
cal simulations is also presented to further support our proposed way to characterize clustering from 
temporal time series at a flxed position. 

PACS numbers: 47.55.Kf, 47.27.T- 
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I. INTRODUCTION 



The phenomenon of preferential concentration of small 
particles and bubbles in turbulent flows attracted much 
attention in recent years, from experimental works |l , [3, 
1l liJEM,0Lto numerical investigations [1, H, [l^, fill . la . 
13 llslJla. [iTt (soj . and theoretical developments [l8| . 
The preferential accumulation is an inertial effect. Par- 
ticles heavier than the fluid are on average ejected from 
vortices, while light buoyant particles tend to accumulate 
in high vorticity regions. Small air bubbles in water (be- 
low 1 mm, typical Reynolds number of order 0(1)) can 
be regarded as a particular kind of non-deformable light 
spherical particles with density negligibly small com- 
pared to the fluid one. In fact, in this size-range, shape 
oscillations or deformations and wake induced effect can 
be reasonably neglected. Strong preferential bubble ac- 
cumulation in core vortex regions is therefore expected 
according to the inertia mechanism, and indeed observed 
experimentally and numerically [l^, [l^] . Next to the 
added mass force and gravity also drag and lift forces 
can affect the clustering. Moreover, the coupling of the 
disperse phase to the fluid flow {two-way coupling) and 
the finite-size effect of particle-particle interaction (4-'way 
coupling) may also result in non-negligible factors of per- 
turbation for preferential concentration of particles and 
bubbles in highly turbulent flows. 

Both the lift force, the two-way coupling, and the four- 
way coupling are notoriously difficult to model in numer- 
ical simulations and a validation of the models against 
numerical simulations is crucial. However, experimen- 
tal measurements on bubbly laden turbulent fiows are 



challenging, as even at very low void fractions 1% in 
volume) the fluid is completely opaque and difficult to ac- 
cess with external optical methods, especially in the bulk 
region of the flow. To experimentally explore the bubble 
clustering in the bulk at high void fraction one there- 
fore has to fall back on intrusive hot-wire anemometer 
measurements. Such measurements had earlier been em- 
ployed to determine the modification of turbulent spectra 
through bubbles [U, [2^ . For the calculation of the 
velocity spectra bubbles hitting the probe had first to 
be identified in the hot-wire signals [2J, [2^ and then fil- 
tered out. In the present paper we employ the very same 
hot-wire time series to obtain information on the bubble 
clustering in the turbulent flow. An alternative method 
to obtain local information on the bubble distribution 
may be phase doppler particle analyzers (26j . 

One could object that measurement from one fixed 
point in space are too intrusive because they can de- 
stroy the clusters, or that they are ineffective in extract- 
ing features of the bubble trapping in turbulent vortical 
structures. The aim of this paper is to demonstrate that 
this is not the case, when using appropriate statistical 
indicators for the analysis of series of bubble colliding 
times on the hot-wire probe. We show that it is pos- 
sible to detect and quantify the micro-bubble clustering 
from a one-point measurement set-up. We compare ex- 
perimental findings with results from numerical simula- 
tions based on Eulerian-Lagrangian approach. Due to 
limitations that we will discuss later, only a qualitative 
agreement among numerics and experiments is expected. 
Nevertheless, we show how this comparison is helpful in 
clarifying the trend in the clustering at changing the tur- 
bulent conditions. 
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II. DETAILS OF THE EXPERIMENT 
METHODS 

The experimental set-up is the Twcntc water channel, 
a vertical duct of square cross section with dimension 
200cTO X 45cTO X 45cm. We refer to Rensen et al. [2l| for 
a detailed description. An array of porous ceramic plates, 
positioned on the top of the channel, is used to generate 
co-flowing small bubbles of average radius, Rq ~ lOOjum, 
as described in [23]. Fluid turbulence is generated by 
means of an active grid, positioned immediately down- 
stream the bubble injection sites. The typical flow is 
characterized by a large mean flow, J7, with turbulent 
fluctuations, u' = {{uz{t) — UY)\^'^, of smaller amplitude. 
The condition u' /U <C 1 assures that Taylor's frozen-flow 
hypothesis can be applied. The dissipative Kolmogorov 
scale measures typically 77 = 400 500 //?ti, while the 
Taylor micro-scale and the integral one, are respectively 
A ~ 30 77, and Lq ~ 500 77. The typical bubble size is of 
the same order, or slightly smaller, than rj. 

We consider microbubble signals extracted from a hot- 
film anemometry probe (300 in size) fixed at the cen- 
ter of the channel. Detection of micro-bubbles is less 
ambiguous than for large bubbles where probe piercing 
and break-up events are highly probable [23|. A micro- 
bubble hitting the probe produces a signal with a clear 
spike. The bubble can be identified by thresholding of 
the velocity time-derivative signal, see Fig. 2 of ref. [23] . 
This identification procedure leads to the definition of a 
minimal cut-off time in the capability to detect clustered 
events, two consecutive bubbles in our records cannot 
have a separation time smaller than r = lG~^sec. Such 
dead-time is mainly linked to the typical response-time 
of the acquisition set-up. Here we consider two time se- 
ries of microbubble measurements, i.e. hitting times, se- 
lected from a larger database because of their uniformity 
and relevant statistics. We will refer to them in the fol- 
lowing as sample (a) and (b). The first sample (a) has 
been taken for a 12 hours long measurement; it consists 
of Ni, ~ 24099 bubbles with a mean hitting frequency 
/ = 0.56 sec^^ . The second sample, (b), is a record of 
11 hours, Nb = 11194 and / ~ 0.28 sec'^ . There are 
two main differences among the experimental conditions 
in which the two samples have been recorded, that is the 
total volume air fraction (called void fraction a), and the 
amplitude of the mean flow and therefore the intensity 
of turbulence. Case (a) has a void fraction of « 0.3% 
and (b) has instead a ~ 0.1%. Note that, even at these 
very small void fractions, the mean number density of 
bubbles amounts to O(IO^) per cubic centimeter. This 
explains the optical opacity of the bulk region of our sys- 
tem. Nevertheless, given the small effect produced by the 
dispersed bubbly phase on the turbulent cascading mech- 
anism [2^ , we consider the discrepancy in a as irrelevant 
for the velocity spectra. In contrast, the difference in the 
forcing amplitude is more important, because it sensibly 
changes all the relevant scales of turbulence as summa- 
rized in Table HI In particular, this leads to different 



values for the minimal experimentally detectable scale: 
lS.rmin — 5?7 for case (a) and Ar^m — 37? for (b), where 
Taylor hypothesis has been used to convert time to space 
measurements, i.e. Ar = t ■ U . In the following, re- 
sults of our analysis will be presented by adopting space 
units made dimensionless by the Kolmogorov scale rj. We 
consider this re-scaling more useful for comparison with 
different experiments and simulations where a mean flow 
may be absent. 



III. DESCRIPTION OF THE STATISTICAL 
TOOLS 

In this section we introduce the statistical tests that 
we will adopt to quantify the clustering. Due to the fact 
that the experimental recording is a temporal series of 
events, we have necessarily to focus on a tool capable to 
identify, from within this one dimensional serie, possible 
signatures of 3D inhomogeneities. 

A first way to assess the presence of preferential con- 
centrations in the experimental records is to compute the 
probability density function (pdf) of the distance, Ar, 
between two consecutive bubbles. Whether the parti- 
cles distribute homogeneously in space, their distribution 
would be a Poissonian distribution and hence the dis- 
tance between two consecutive bubbles would be given 
by the well know exponential expression: /9exp(— pA7'); 
where p ~ f /U is the number of bubbles per unit length 
(i.e. their density) [l^. Due to the presence of turbu- 
lence we expect that in general the spatial distribution of 
the bubbles will differ from a Poissonian distribution: in 
any case it is natural to expect that for separation-scales 
large enough the exponential form of the pdf should be 
recovered. In fact, pairs of successive bubbles with large 
separations Ar, larger then any structures in the flow, 
are expected to be uncorrelated, memory-less, events. 

Due to the possible accumulation on small scales (clus- 
tering of bubbles) the long tail of the pdf may have an 
exponential decay rate that is different from the global 
mean, p. The tail of the experimentally measured pdf 
can be fitted with an exponentially decaying function, 
A ■ exp(— p/jAr), with a rate that we call ph, where h 
stands for homogeneous. In case of small-scale clustering 
we expect ph to be smaller than p. As an indicator of the 
fraction of bubbles accumulated in turbulent structures, 
we use the coefficient C = 1 — p^/ p, whose value varies 
in between and 1. 

The test so far introduced is useful but only provides an 
indication on how homogeneously distributed the bubbles 
are at small scales, while it gives no indication on their 
possible "large-scale" correlations. Here we introduce a 
second, more comprehensive, statistical test particularly 
convenient to reveal the scales at which the inhomogene- 
ity develops. The idea is to compute the coarse-grained 
central moments of the number of bubbles, on a window 
of variable length r, pP = {{n — {n)rY)r- The length 
of the window r will be the scale at which we study 
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Lo (cm) 


U (cm/s) 


u (cm/s) 


Rex 


Teddy (s) 


Tri (ms) 




Ujj (mm/s) 


Ret 


Ro/ri 


St 


gn/ur, 


a) 


22.6 


19.4 


1.88 


206 


12.0 


151. 


388. 


2.57 


4.4 


0.26 


0.007 


4.2 


(b) 


23.1 


14.2 


1.39 


180 


16.6 


240. 


489. 


2.04 


4.4 


0.20 


0.004 


5.3 



TABLE I: Relevant turbulent scales and bubble characteristics for the two experimental samples analyzed. Fluid turbulent 
quantities have been estimated from one-dimensional energy spectra. From left to right: integral scale, Lo, mean velocity, U, 
single-component root mean square velocity, u' , Taylor Reynolds number. Rex, large eddy turnover time, Teddy, dissipative 
time (r,,) space (ry) and velocity (u,,) scales, bubble Reynolds number, Reb, bubble-radius and Kolmogorov-length ratio, Ro/r], 
Stokes number, St, ratio between terminal velocity in still fluid and dissipative velocity scale, gTt/un- 



whether the distribution resembles a homogeneous one. 
We will focus on scale dependent kurtosis and skewness 
excesses, respectively: K{r) = — 3 and S{r) = 

fj,f/{fi^)^^^. A random distribution of particles spatially 
homogeneous with mean density p corresponds to the 
Poissonian distribution : p{n) = exp(— p r)(pr)"(n!)^^, 
where r is the length of the spatial window and n is the 
number of expected events. Therefore, once the par- 
ticle space rate p is given, the value of any statistical 
moment can be derived for the corresponding window 
length r. A spatially Poissonian distribution of parti- 
cles implies the functional dependences K{r) — (pr)^^ 
and S{r) = {pr)~^^^ . Furthermore, we note that at the 
smallest scale, when r = Armin, we reach the singular 
limit (shot-noise limit) where for any given space- window 
we can find none or only one bubble and all statistical mo- 
ments collapse to the same value. This latter limit, which 
is by the way coincident with Poisson statistics, repre- 
sents our minimal detectable scale. We are interested 
in departures from the shot-noise/random- homogeneous 
behavior for the statistical observables K{r) and S{r). 



IV. RESULTS OF THE ANALYSIS ON 
EXPERIMENTAL DATA 

In figure [T] we show the computed pdf(Ar) for the 
two data samples considered. Deviations from global 
homogeneity are clear if the shape of the histogram 
is compared with the solid line representing the pdf 
/9exp(— pAr). These deviations are slightly more pro- 
nounced in the more turbulent case (a) as compared to 
case (b). Nevertheless, one can notice that the pure ex- 
ponentially decaying behavior, i.e. homogeneity, is re- 
covered from distances of the order of ©(lOOyy) up to the 
large scales. The dotted line on Fig. [U which represents 
the linear fit on the long homogeneous tail in the inter- 
val [10^,2 ■ 10'^]?7, and the inset boxes, where the pdf is 
compensated by the fit, shows this latter feature. The 
evaluation of the coefhcient C leads to values for the rela- 
tive bubbles excess in clusters corresponding to 19 % for 
case (a) {Rex ^ 206) and 10 % for case (b) (Re), ~ 180), 
confirming the trend of stronger concentration in flows 
with stronger turbulence level. In figure [5] we show the 
Kurtosis and Skewness behavior, evaluated for the two 
cases (a)-(b), in a comparison with the Poissonian de- 
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FIG. 1: Probability density function of distance between 
successive bubbles, pdf(Ar). Exponential behavior, pe~''^^, 
(solid line) and exponential fit, A ■ e~'''^^^ , of the large-scale 
tail (dashed line) are reported. The inset shows the pdf(Ar) 
compensated by the fitted large-scale exponential behavior, 
i.e. the pdf(Ar) divided by A ■ e'"'-^'' . 



pendence. We observe in both cases a clear departure 
at small scale from the scaling implied by the global 
homogeneity, which is only recovered at the large scale 
(> Lo — 5OO77) where the data points falls roughly paral- 
lel to the Poisson line. The departure from the Poisson 
line, that is noticeable already at the scales immediately 
above Ar^in , is an indication that bubbles form clusters 
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FIG. 2: Scale dependent Kurtosis, K(r), for case (a) (top) and 
(b) (bottom). Dotted lines represent the Poissonian behavior, 
that is i4'(p)(r) = (p r)~^. Notice that the Poisson scaling 
behavior is reached for large r-windows only scaling wise. In 
the insets the scale dependent Skewness, S{r), behavior is 
shown. Again the Poissonian relation is drawn S(^p){r) — 
[p r)^^''^ (dotted line). 



even at the smallest scale we are able to detect, that is 
even below 5r/ for case (a) or 877 for case (b) . We observe 
that for the less turbulent case, (&), the departure from 
the homogeneous scaling is less marked. A comparison 
with synthetic Poisson samples of an equivalent number 
of bubbles, that we have tried, shows that the available 
statistics is sufficient to validate the deviations from the 
homogeneity discussed so far. Scale dependent deviation 
from Poisson distribution is an evidence of the fact that 
the dispersed microbubblcs are trapped within the dy- 
namical vortical structures of turbulence. Furthermore, 
we observe that gravity plays a minor role in this dynam- 
ics. In fact, on average the bubbles are swept down by 
the mean flow and gTb/ujj ~ 0{1) (see Tab. |l|, which 
implies that even the smallest vortical structures of the 
flow may trap bubbles [l^. Therefore, it is mainly the 
inertia that drives the bubble accumulation in the flow. 



V. RESULTS OF THE ANALYSIS ON 
NUMERICAL DATA 

To give further evidence for the robustness of the sug- 
gested statistical analysis of the hot-wire time series, we 
now repeat the very same procedure with numerical simu- 
lation data. We employ standard numerical tools already 
described and discussed in details in In short, we 

integrate Lagrangian point-wise bubbles evolving on the 
background of an Eulcrian turbulent fleld. The equation 
for the evolution of the point-wise bubble is the following: 

dv Du 1 , , , , , , 

— = 3-7;^ {v-u)-2g- {v~u)xu} (1) 

at at Tfc 

where u and u) are respectively the fluid velocity and vor- 
ticity computed at the bubble position and constitute a 
simplified version of the model suggested in [2^ . The Eu- 
lerian flows is a turbulent homogenous and isotropic fleld 
integrated in a periodic box, of resolution 128'^, seeded 
approximately with 10^ bubbles, corresponding to a void 
fraction a = 4.5%. Since previous experimental and 
numerical studies [l^, [l^l have revealed that the effect 
of bubbles on strong unbounded turbulence is relatively 
weak, our numerical bubbles are only coupled in one- 
way mode to the fluid, i.e. bubbles do not affect the 
fluid phase. The bubble-Reynolds number Ret is set to 
unity and the Stokes number is St ^ 1. Therefore, the 
bubble radius is of order 77, and the bubble terminal ve- 
locity vt = "^gTb in still fluid is smaller than the smallest 
velocity scale u^. As (and actually even more than) in 
the experiment, the role of gravity is marginal. In Table 
im we report details of the numerical simulations, these 
are chosen trying to match the experimental numbers. 
However, we could not reach the same scale separations 
as in the experiments. In the bottom panel of TabHIl we 
translated the numerical units to their physical equiva- 
lent. We note that in the numerics the Stokes number, 
St = Tb/Tjj, which is an indicator of the degree of bubble 
interaction with turbulence, can not be as low as in the 
experiments. To achieve the same St would require too 
much CPU time. For practical reason the Stokes values 
adopted in our numerics are roughly one order of mag- 
nitude larger than in the experiments, although always 
much smaller than unity, St 1. Under this conditions, 
simple spatial visualization [l^ show strong bubble accu- 
mulation in nearly one-dimensional elongated structures 
in correspondence to high enstrophy regions (identifled 
as vortex fllaments. As already stated, our goal is to use 
the numerics to confirm the behavior of the suggested 
observables. To this end, we put 128 virtual pointwise 
probes in the flow and recorded the hitting times of the 
bubbles, which we give a virtual radius Rq. The bubble 

radius is related to the bubble response time rt,, namely 
1/2 

Tf, = i?o = (STfei') when assuming no-slip boundary 
conditions at the gas-liquid interface. 

An important difference between the experiments and 
the numerics is the mean flow: It is present in the experi- 
ment while intrinsically suppressed in the simulations. In 
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St 




(a') 


5.0 


1.4 


94 


3.6 


0.093 


0.025 


0.275 


1.0 


1.13 


0.14 


0.55 


(b') 


5.0 


1.0 


87 


4.9 


0.147 


0.032 


0.218 


1.0 


0.89 


0.09 


0.69 




Lo (cm) 


u (cm/s) 


Rex 


Teddy (ms) 


Tr, {ms) 


V {^J■'m) 


Un (cm/s) 


Reb 


Ro/v 


St 


gn/u,, 


(a') 


0.41 


7.2 


94 


57.3 


4.7 


68.7 


1.45 


1.0 


1.13 


0.14 


0.55 


(b') 


0.46 


5.5 


87 


82.5 


7.3 


85.7 


1.16 


1.0 


0.89 


0.09 


0.69 



TABLE II: Relevant turbulent scales and bubble characteristics for the two numerical simulation performed. The top part 
reports the actual values in numerical units from the simulation, the bottom part shows for comparison the corresponding 
physical equivalent quantities for air bubbles in water, this is to better appreciate similarities/differences with the experimental 
conditions of Table |l] The values on the bottom part are computed starting from the dimensionless quantities Rex, Reb, St, 
and by assuming i/ = 10~® m^ and g = 9.8 m s^^. 
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FIG. 3: Numerical result on the probability density function 
of distance between successive bubbles, pdf(Ar). Case (a') 
(top) is the most turbulent. In the inset the same compen- 
sated plot as in Fig. [T] 




FIG. 4; Numerical result on scale dependent kurtosis, K{r), 
for case (a') (top) and (b') (bottom), and Poissonian behavior 
(dotted). In the insets the scale dependent skewness, S{r), 
behavior is shown. 



the numerical simulations the time is connected to space 
displacements through the relation Ai? = At • it', where 
u' is the root mean square velocity. The level of turbu- 
lence, given the available resolution, has been pushed as 
high as possible {Re\ ~ 90) to obtain a better analogy 
with the experiment. Also in the numerical simulations 



two cases with different Reynolds numbers are consid- 
ered, see again Table HIl 

In Figures [3] and [5] we show the results of the statistical 
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analysis of clustering from the time series obtained from 
the numerical virtual probes. These two figures should 
be compared with the analogous experimental findings 
already discussed and shown in Figures [T] and [21 Some 
qualitative similarities are striking. First, starting from 
Fig. [3l we observe that deviations from random and ho- 
mogeneous, i.e. pure exponential behavior, are relevant 
at small scales. This feature is confirmed by the scale 
dependent kurtosis and skewness of Fig. [4l where depar- 
ture from the Poisson scaling already starts below 77 scale. 
Second, the most turbulent case is the most clusterized, 
(a') {Rex ~ 94) more than (b') {Rex 87). The evalu- 
ation of the fraction of clustered bubbles, based on the 
fit of the pdf(Ar) as in the experiment, gives the value 
29% for (a') and 37% for (b'). Though the quahtative be- 
havior of the statistical indicators is the same, also some 
important differences arise in this comparison. First of 
all, full homogeneity in the numerics seems to be recov- 
ered already at scales of order O{10ri), whereas in the 
experiments if was only recovered at ©(lOOjy). Further- 
more, the deviations from Poisson distribution and the 
fraction of clustered bubbles are definitely stronger in the 
numerics. There are several possible interpretation for 
this mismatch, including the possible incompleteness of 
the employed model eq. First, some physical effects 
have been neglected: the fluid-bubble and the bubble- 
bubble couplings and the associated finite size effects (in 
the present conditions bubbles can overlap!). A second 
reason can be the different degree of bubble interaction 
with turbulence, a quantity that is parametrized by the 
Stokes number St = Tt,/T,i. The estimated St in the ex- 
periment is roughly one order of magnitude smaller than 
in the simulation. This corresponds to bubbles that react 
faster to the fluid velocity changes and hence to bubbles 
that closely follow the fluid particles and accumulate less. 
Such a trend is also confirmed by our numerics. 



erential bubble concentration from single-point hot-wire 
anemometer measurements in the bulk of a turbulent 
channel. Our tools clearly show that the experimental 
records display bubble clustering. The fraction of bub- 
bles trapped in such structures is indeed considerable 
and can be estimated to be of the order of 10%. The 
scale-dependent deviations from random homogeneous 
distribution, that we associate to typical cluster dimen- 
sion, extends from the smallest detectable scale, 0{ri), 
to scales in the low inertial range, O{100r]). Accumula- 
tion of bubbles is enhanced by increasing the turbulence 
intensity. Comparison with present Eulerian-Lagrangian 
simulations, where point-like bubbles strongly accumu- 
late in vortex core regions, shows similar qualitative fea- 
tures and trends. 

We hope that our explorative investigation will stim- 
ulate new dedicated experiments and numerical simula- 
tions to further quantify the clustering dynamics as func- 
tion of Reynolds number and particle size, type, and 
concentration. The challenge is to further develop and 
employ quantitative statistical tools to allow for a mean- 
ingful comparison between experiment and simulations, 
in order to validate the modeling of particles and bubbles 
in turbulent flow. 
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VI. CONCLUSIONS 



We have performed statistical tests in order to detect 
and quantitatively characterize the phenomenon of pref- 
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